function [likel,ssvec,flag_ok,lht]=likelCumSplit(parest,parvec,parposest,funcmod,Y,trainvec,solveopt,addsol)
% =========================================================================
%% LikelNanSplit.m
%
% This code estimates the likelihood of 
% 1) NaN model 
% 2) Mixed Frequency 
% 3) One time Split in Solution matrices, where all matrices (GG,RR,C,SDX,ZZ) are allowed to change in the second part of
% the sample at a time T*
%
%% *Additional inputs needed*
% *addsol.dateBegSecond*  Numeric Calendar date when second sample begins
% (provided in settings)
% *addsol.rowBegSecond*  Position in the sample vector where second sample
% begins (computed in LMJ_initial.m)
%
% Matrices of the second sample are given in the structure second, usually
% called FLAG, in the solution of the model, e.g.
%
% [GG, RR, CONS, eu, SDX, ZZ, initss, ssvec, flag, ssnames]=feval(funcmod,param,options,addsol);
%
% =========================================================================
%% Initialize parameters
likel=-1e40; flag_ok=0;
parvec(parposest)=parest;
%% Model solution, second sample stored in structure second (second sample) 
[G,R,C,eu,SDX,Z,structOne,ssvec,structTwo]=feval(funcmod,parvec,solveopt,addsol);

% A. Non-existence in either sample
if isequal(eu,[1;1])==0 || field.euok==0 || isequal(structOne.eu,[1;1])==0
    return
end

% =========================================================================
%% 2. Set-up extended states space 
%     See buildState.m by AB for an alternative 

%% 2.1 Dimensions and time invariant matrices 
[T,ny] = size(Y); 
ns      = size(G,1);
nx   = size(SDX,2);

if size(Z,1)~=structOne.nzHF;error('Rows in Z must match number of variables observed in High frequency'); end 
if size(structOne.Z,1)~=structOne.nzHF;error('Rows in Z2 must match number of variables observed in High frequency'); end 


Z1tilda  = [Z zeros(structOne.nzHF,structOne.nzLF); ...
    zeros(structOne.nzLF,ns) eye(structOne.nzLF)];
Z2tilda =  [structOne.Z zeros(structOne.nzHF,structOne.nzLF);...
    zeros(structOne.nzLF,ns) eye(structOne.nzLF)];

% Constrant vector 
C1   = [C; structOne.A*C];
C2 = [structOne.C; structOne.A*structOne.C];

% casevec: vector determining the particular page of the T and R time 
% varying matrices
casevec =kron( ones(T/structOne.NUMIT,1) , (1:structOne.NUMIT)' );

% type of temporal aggregation
if structOne.flag_tagsum == 1
    sigma =[0 ones(1,structOne.NUMIT-1)];
else
    sigma =(1:structOne.NUMIT);
end


%% 2.2 Page matrices with pages corresponding to quarters {1,2,3,4} 
Tpag   =  zeros(ns+structOne.nzLF,ns+structOne.nzLF,structOne.NUMIT);
Tpag2  =  zeros(ns+structOne.nzLF,ns+structOne.nzLF,structOne.NUMIT);
Rpag   =  zeros(ns+structOne.nzLF,  nx    ,structOne.NUMIT);
Rpag2  =  zeros(ns+structOne.nzLF, nx    ,structOne.NUMIT);
for ii=1:structOne.NUMIT 
    if structOne.flag_tagsum == 1
        Tpag(:,:,ii) = [G zeros(ns,structOne.nzLF);...
            structOne.A*G eye(structOne.nzLF)*sigma(ii)];
        Rpag(:,:,ii) = [R; structOne.A*R];
        Tpag2(:,:,ii) = [structOne.G zeros(ns,structOne.nzLF);...
            structOne.A*structOne.G eye(structOne.nzLF)*sigma(ii)];
        Rpag2(:,:,ii) = [structOne.R; structOne.A*structOne.R];
    else
        Tpag(:,:,ii) = [G zeros(ns,structOne.nzLF); ...
            (1/sigma(ii))*structOne.A*G eye(structOne.nzLF)*(sigma(ii)-1)/sigma(ii)];
        Rpag(:,:,ii) = [R; (1/sigma(ii))*structOne.A*R];
        Tpag2(:,:,ii) = [structOne.G zeros(ns,structOne.nzLF);...
            (1/sigma(ii))*structOne.A*structOne.G ...
            eye(structOne.nzLF)*(sigma(ii)-1)/sigma(ii)];
        Rpag2(:,:,ii) = [structOne.R; ...
            (1/sigma(ii))*structOne.A*structOne.R];
    end 
end 


%% Demeaning using the split sample
%  First Sample (C1)
if any(C1~=0)==1
    Y(1:addsol.rowBegSecond-1,:)=...
        Y(1:addsol.rowBegSecond-1,:)-repmat((Z1tilda*C1)',[addsol.rowBegSecond-1 1]);
end 
% Second sample (second.second.C)
if any(C2~=0)==1;
    Y(addsol.rowBegSecond:end,:)=...
        Y(addsol.rowBegSecond:end,:)-repmat((Z2tilda*C2)',[(T-addsol.rowBegSecond+1) 1]);
end 
Y=Y';

%% Initialization 
pshat=feval(@lyapunov_symm,Tpag(:,:,casevec(1)),...
    Rpag(:,:,casevec(1))*(SDX')*(Rpag(:,:,casevec(1))*(SDX')'));
shat=zeros(ns+structOne.nzLF,1);

%% Forward Filter

% storage matrices
lht =zeros(T,1);
Zdim=zeros(T,1); 
W   =eye(ny);

% Filter, first part
for ii=1:addsol.rowBegSecond-1;
    
    % Handling of missing observations 
    ytt=Y(:,ii);
    
    % Determine W and position of the NAN  
    ind =~isnan(ytt);
    
    ytt=ytt(ind);    
    Zdim(ii)=length(ytt); 
    Ztt=W((ind==1),:)*Z1tilda;
    [shat,pshat,lht(ii)]=...
        feval(@kf,ytt,Ztt,shat,pshat,Tpag(:,:,casevec(ii)),Rpag(:,:,casevec(ii))*(SDX'));
end

% Filter, second part
for ii=addsol.rowBegSecond:T;
    
    % Handling of missing observations 
    ytt=Y(:,ii);
    
    % Determine W and position of the NAN  
    ind =~isnan(ytt);
    
    ytt=ytt(ind);    
    Zdim(ii)=length(ytt); 
    Ztt=W((ind==1),:)*Z2tilda;
    [shat,pshat,lht(ii)]=...
        feval(@kf,ytt,Ztt,shat,pshat,Tpag2(:,:,casevec(ii)),Rpag2(:,:,casevec(ii))*(structOne.SDX'));
end

%% Likelihood plus integration constant

lht=lht(trainvec(1):trainvec(2));
likel=-0.5*log(2*pi)*(sum(Zdim))+sum(lht);
flag_ok=1;

%% End of file
end